
import random

###-------------------------------------###
###### Polynomial Multiplication  ######
###-------------------------------------###
import collections
def mulPol(adict, bdict):
    """
    >>> a, b = [{0:1, 1:1}] * 2
    >>> mulPol(a, b)
    Counter({1: 2, 0: 1, 2: 1})
    >>> sorted(reduce(mulPol, [{0:1, 1:1}]*3, {0:1}).items())
    [(0, 1), (1, 3), (2, 3), (3, 1)]
    """
    ans = collections.Counter()
    for k,v in adict.items():
        for k1,v1 in bdict.items():
            ans[k+k1] += v * v1
    return ans

import lxc
def addPol(adict, bdict):
    """
    >>> a, b = [{0:1, 1:1}] * 2; sorted(addPol(a, b).items()) == [(0, 2), (1, 2)]
    True
    >>> b={0:1,1:1,2:3}; sorted(addPol(a, b).items()) == [(0, 2), (1, 2), (2, 3)]
    True
    """
    #>>> b[2]=3; sorted(addPol(a, b).items()) == [(0, 2), (1, 2), (2, 3)]
    import lxc
    return lxc.dtUnion(adict, bdict, lambda a,b: a+b)

class ExpoPol(dict):
    """
    This class is to support the Exponential Generation Functions.
    """
    def __init__(self, *argv, **argvv):
        super(ExpoPol, self).__init__()
        self.update(*argv, **argvv)
    def __repr__(self):
        #return "ExpoPol(%s)" % (dict.__repr__(self))
        return "ExpoPol(%s)" % (super(ExpoPol, self).__repr__())
    def __str__(self):
        a = sorted(self.items())
        return " + ".join(["%fx^%d/%d!"%(j*fac(i), i, i) for (i, j) in a])
    def sstr(self):
        a = sorted(self.items())
        return " + ".join(["%dx^%d/%d!"%(j*fac(i), i, i) for (i, j) in a])
    def getPowerItem(self, k):
        return self[k] * fac(k)
    @staticmethod
    def create(string):
        d = Pol.create(string)
        for i in d.keys(): d[i] = d[i] / float(fac(i))
        return ExpoPol(d)

class Pol(dict):
    """
    This class is to support the Normal Generation Functions.
    """
    def __str__(adict):
        a = sorted(adict.items())
        return " + ".join(["%dx^%d"%(j, i) for (i, j) in a])
    def __repr__(self):
        return "Pol(%s)" % (super(Pol, self).__repr__())
    @staticmethod
    def analyse(string):
        tmp = "".join([i for i in string if i in 'x+-0123456789^'])
        x = tmp.split("+")
        ans = collections.Counter()
        for (v, k) in [map(int, i.split("x^")) for i in x]: ans[k] += v
        return ans
    @staticmethod
    def create(string):
        tmp = "".join([i for i in string if i in 'x+-0123456789'])
        x = tmp.split("+")
        ans = collections.Counter()
        for (v, k) in [map(int, i.split("x")) for i in x]: ans[k] += v
        return Pol(ans)
    

###-------------------------------------###
###### Combin. Math Impl. and Test   ######
###-------------------------------------###

def A(n, r):
    """
    >>> A(5, 3)
    60
    """
    ans = 1
    for i in xrange(n-r+1, n+1): ans *= i
    return ans

def C(n, r):
    """
    >>> C(5, 2) == C(5, 3)
    True
    """
    ans = 1
    for i in xrange(n-r+1, n+1): ans *= i
    for i in xrange(1, r+1): ans /= i
    return ans

def fac(n):
    """
    fac(5) = 5 * 4 * 3 * 2 * 1 = 120
    >>> all([fac(5)==120, fac(0) == 1])
    True
    """
    ans = 1
    for i in xrange(2, n+1): ans *= i
    return ans

def fib(n, init0=1, init1=1):
    """
    get fib n-term
    >>> map(fib, range(9))
    [1, 1, 1, 2, 3, 5, 8, 13, 21]
    """
    for i in xrange(n-2): init0, init1 = init1, init0 + init1
    return init1

import math
### Todo: How to make it simple?
def ln(n):    
    return math.log(n)

###-------------------------------------###
###### Statis. Util. Impl. and Test  ######
###-------------------------------------###

class Acc():
    """
    >>> stda.addl(range(10)), stda.gavg() == 4.5, stda.reset()
    (None, True, None)
    """
    def __init__(self, initList = None):
        self.reset(initList)
    def reset(self, alist=None):
        if alist == None:
            self.alist = list()
        else: self.alist = alist
    def add(self, num):
        self.alist.append(num)
    def __lshift__(self, num):
        self.add(num)
    def addl(self, alist):
        for i in alist: self.add(i)
    def dlt(self, num):
        self.alist.remove(num)
    def dlt(self):
        self.alist.pop()
    def gsum(self):
        return sum(self.alist)
    def gavg(self):
        return (sum(self.alist)+0.0) / len(self.alist)

stda = Acc()

########################################################
########################################################
########################################################

import fractions
def gcd(a, b):
    """
    >>> gcd(12, 18)
    6
    """
    return fractions.gcd(a, b)

def lcm(a, b):
    """
    >>> lcm(12, 18)
    36
    """
    return (a*b) / gcd(a, b)


########################################################
########################################################
########################################################

def power(a, b):
    """
    Actually, the performance is worse.
    The native ** implementation is no doubt better.
    It is only for an algorithm illustration.
    """
    b = int(b)
    s, t = 1, a
    while (b > 0):
        if (b & 1L == 1):
            s *= t
        t *= t
        b >>= 1
    return s

def powerMod(a,b,c):
    """
    It is just for the scenaria when b is positive integer
    >>> native = lambda a, b, c: a ** b % c
    >>> all(map(lambda a:native(*a)==powerMod(*a), [(3,21,38), (4, 19, 97)]))
    True
    """
    b = int(b)
    s, t = 1, a
    while (b > 0):
        if (b & 1L == 1):
            s = (s * t) % c
        t = (t*t) % c
        b >>= 1
    return s
    
########################################################
########################################################
########################################################

### Number Theory

#Very slow
def getPrimeList(num):
    '''
    >>> getPrimeList(25), getPrimeList(50) == [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
    ([2, 3, 5, 7, 11, 13, 17, 19, 23], True)
    '''
    alist = range(3, num+1, 2)
    sqrt = int(math.sqrt(num))
    k = 0
    while (k < len(alist)):
        e = alist[k]
        if e > sqrt: break
        alist = alist[0:k+1] + filter(lambda a:a % e != 0, alist[k+1:])
        k += 1
    return [2] + alist

import math

#Fastest without numpy
def primeGen(num):
    '''
    >>> map(lambda a: primeGen(a) == getPrimeList(a), [25, 50])
    [True, True]
    '''
    f = [1] * (num+1)
    sqrt = int(math.sqrt(num))
    for s in xrange(3, sqrt+1, 2):
        if f[s] == 1:
            for multi in xrange(s+s, num+1, s):
                f[multi] = 0

    res = [(i<<1)+1 for i, v in enumerate(f[1::2], 0) if v == 1]
    res[0] = 2
    return res

###
def iprimeGen(num):
    '''
    >>> map(lambda a: list(iprimeGen(a)) == getPrimeList(a), [25, 50])
    [True, True]
    '''
    yield 2
    f = [1] * (num+1)
    sqrt = int(math.sqrt(num))
    for s in xrange(3, sqrt+1, 2):
        if f[s] == 1:
            yield s
            for multi in xrange(s+s, num+1, s):
                f[multi] = 0
    s += 2
    for i, v in enumerate(f[s::2], s/2):
        if v == 1: yield ((i<<1)+1)

#Fastest of all
def getNpPrimeList(num):
    '''
    sumAnd(map(lambda a: getNpPrimeList(a) == getPrimeList(a), [25, 50, 100]))
    True
    '''
    num = num+1
    a = np.arange(num)
    a[0:2] = -1
    a[4:num:2] = -1
    s = 3
    sqrt = int(math.sqrt(num))
    while s <= sqrt:     
        if (a[s] != -1):
            a[s+s:num:s] = -1
        s += 2
    return np.compress(np.greater(a, -1), a).tolist()

def gPrimeL(num):
    '''
    >>> sumAnd(map(lambda a: gPrimeL(a) == getPrimeList(a), [25, 50, 100]))
    True
    '''
    try:
        import numpy as np
        return getNpPrimeList(num)
    except:
        return primeGen(num)

def isPrimeRisk(num):
    a = [2, 3, 5, 7, 11, 13, 17, 19, 23]
    if num<=23: return num in a
    a = [random.randint(2, num-1) for i in xrange(50)]
    for e in a:
        if powerMod(e, num-1, num) != 1:
            return False
    return True

#Very slow
def getPrimeListRisk(num):
    '''
    >>> getPrimeListRisk(25)
    [2, 3, 5, 7, 11, 13, 17, 19, 23]
    >>> getPrimeListRisk(100)
    [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
    >>> getPrimeListRisk(101)
    [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101]
    '''
    ans = []
    for i in xrange(2, num+1):
        if isPrimeRisk(i):
            ans.append(i)
    return ans

# Work for num > 1
def isPrime(num):
    '''
    >>> primelist = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
    >>> sigor << map(isPrime, filter(lambda e: e not in primelist, range(2, 50)))
    False
    >>> sumAnd(map(isPrime, primelist))
    True
    '''
    for i in xrange(2, int(math.sqrt(num))+1):
        if (num % i == 0):
            return False
    return True

### Implementation is NOT efficient. Try Time Analysis
def getFactors(num):
    a = set()
    for i in xrange(2, int(math.sqrt(num))+1):
        if (num % i == 0):
            a.add(i)
            a.add(num/i)
    return sorted(a)

def getPrimeFactors(num):
    if num == 1: return []
    if num == 2: return [2]
    c = getPrimeList(num+1)
    a = []
    k = 0
    while (num != 1):
        if num % c[k] == 0:
            num /= c[k]
            a.append(c[k])
        else: k+=1
        if (k >= len(c)):
            break
    return a


class DecimalComparator():
    """
    >>> [1.8 + 0.1 == 1.9, stdf.equ(1.8+0.1, 1.9)]
    [False, True]
    """
    def __init__(self, eps = 10e-8):
        self.eps = eps
    def g(self, a, b):
        return a - b > self.eps
    def ge(self, a, b):
        return a > b - self.eps
    def l(self, a, b):
        return a < b - self.eps
    def le(self, a, b):
        return a < b + self.eps
    def equ(self, a, b):
        return abs(a-b) < self.eps

stdf = DecimalComparator()

### sum --- Function Extention

import operator
def product(alist):
    return reduce(operator.mul, alist, 1)

def sumAnd(alist):
    return reduce(operator.and_, alist, True)

def none(alist):
  #  return all(map(lambda a:not a, alist))
  #  return all(map(operator.not_, alist))
  #  return not(reduce(operator.or_, alist, False))
    '''
    >>> alist = [False] * 10; blist = [False]*10 + [True] + [False]
    >>> none(alist) == all(map(operator.not_, alist)), none(blist) == (not (reduce(operator.or_, blist, False)))
    (True, True)
    >>> naive = lambda alist:all(map(operator.not_, alist))
    >>> all(map(lambda a:naive(a)==naive(a), [alist, blist] ))
    True
    '''
    for i in alist:
        if i == True:
            return False
    return True
### Same
### sigand = lshiftFun1(sumAnd)
class SUM_AND():
    def __lshift__(self, alist):
        return sumAnd(alist)
sigand = SUM_AND()

def lshiftFun1(fun):
    class EXAMPLE():
        def __lshift__(self, e):
            return fun(e)
    return EXAMPLE()

def sumOr(alist):
    return reduce(operator.or_, alist, False)
sigor = lshiftFun1(sumOr)

def sumXor(alist):
    return reduce(operator.xor, alist, False)
sigxor = lshiftFun1(sumXor)

### Median Function

def _median(alist, percentage=0.5):
    nth = len(alist) * percentage
    kth = int(nth)
    if (nth-kth<1e-18):
        return (alist[kth] + alsit[kth+1]) / 2, kth
    else: return alist[kth+1], kth

def median_in(alist, percentage=0.5):
    alist.sort()
    return _media(alsit, percentage)
    
def median(alist, percentage=0.5):
    return (sorted(alist), percentage)

def _ridlist0(alist, percentage=0.5):
    nth = int(len(alist) * percentage / 2)
    return nth, alist[nth:len(alist)-nth]


def ridlist0_in(alist, percentage=0.5):
    alist.sort()
    return _ridlist0(alist, percentage)

def ridlist0(alist, percentage=0.5):
    return ridlist0(sorted(alist), percentage)

def _ridlist1(alist, percentage=0.5):
    nth = int(len(alist) * percentage / 2)
    return alist[:nth], alist[nth:len(alist)-nth], alist[len(alist)-nth:]

def ridlist1_in(alist, percentage=0.5):
    alist.sort()
    return _ridlist1_in(alist, percentage)

def ridlist1(alist, percentage=0.5):
    return _ridlist1(sorted(alist), percentage)

### ###
# Vector Section
### ###

sqr = lambda e: e*e

class Vector:

    @classmethod
    def norm(self, a): return math.sqrt(self.norm2(a))
    @classmethod
    def calR(self, a, b): return math.sqrt(self.calRR(a, b))
    @classmethod
    def center(self, alist):
        """
        >>> Vector2.center([(2, 3), (4, 5), (3, 7)]) == (3, 5)
        True
        """
        t = self.getZero()
        for e in alist:
            t = self.add(t, e)
        t = self.mul(t, 1.0 / len(alist))
        return t
    @classmethod
    def relative(self, e, center):
        return self.sub(e, center)
    @classmethod
    def centralize(self, loc, center=None):
        if center == None:
            center = self.center(self.loc)
        return [self.sub(l, center) for l in loc], center
    # Caution: That's old center, New center is Zero
    # this only used to measure the movement vector

    @classmethod
    def centralize_scale(self, loc, ratio=1.2, center=None):
        """
        >>> v = Vector2
        >>> v.centralize_scale([(6, 0), (-4, 0), (1, 5), (1, -5), (1, 0)])==[(6, 0), (-6, 0), (0, 6), (0, -6), (0, 0)]
        True
        >>> v.centralize_scale([(6, 1), (-4, 1), (1, 6), (1, -4), (1, 1)])==[(6, 0), (-6, 0), (0, 6), (0, -6), (0, 0)]
        True
        """
        if center == None: center = self.center(loc)
        return [self.mul(self.sub(l, center), ratio) for l in loc]

    @classmethod
    def scale(self, loc, ratio=1.2, center=None):
        """
        >>> v = Vector2
        >>> v.scale([(5, 0), (-5, 0), (0, 5), (0, -5)])==[(6, 0), (-6, 0), (0, 6), (0, -6)]
        True
        >>> v.scale([(6, 0), (-4, 0), (1, 5), (1, -5), (1, 0)])==[(7, 0), (-5, 0), (1, 6), (1, -6), (1, 0)]
        True
        >>> v.scale([(6, 1), (-4, 1), (1, 6), (1, -4), (1, 1)])==[(7, 1), (-5, 1), (1, 7), (1, -5), (1, 1)]
        True
        """
        if center == None: center = self.center(loc)
        return [self.add(self.mul(self.sub(l, center), ratio), center) for l in loc]
        #ans = []
        #for l in loc:
        #    v = self.sub(l, center)
        #    newv = self.mul(v, ratio)
        #    ans.append(self.add(newv, center))
        #return ans
        #OR
        # if center == None: center = self.center(loc)
        # return map(self.cenralize_scale(loc, ratio, center), lambda e:self.add(e, center))
    @classmethod
    def toR(self, vec, tgtr):
        ratio = tgtr / self.norm(vec)
        return self.mul(vec, ratio)
    @classmethod
    def floorR(self, vec, floorR):
        """
        >>> v = Vector2
        >>> v.floorR((6, 8), 15) == ((9, 12))
        True
        """
        r = self.norm(vec)
        if r<floorR: return self.mul(vec, floorR / r)
        else: return vec
    @classmethod
    def ceilR(self, vec, ceilR):
        """
        >>> v = Vector2
        >>> v.ceilR((6, 8), 5) == (3, 4)
        True
        """
        r = self.norm(vec)
        if r>ceilR: return self.mul(vec, ceilR / r)
        else: return vec    
    
class Vector2(Vector):
    """
    >>> v = Vector2
    >>> a, b = (3, 4), (0, 0); Vector2.norm2(Vector2.sub(a, b)) == 25
    True
    >>> v.add(a, b) == a
    True
    >>> v.mul(a, 3) == (9, 12)
    True
    >>> v.dot(a, a) == (9, 16)
    True
    >>> v.neg(a) == (-3, -4)
    True
    >>> v.dot(a, b) == (0, 0)
    True
    >>> lll = lambda e0, e1: v.norm2(v.sub(e0, e1)) == v.calRR(e0, e1)
    >>> test = [[3, 1, 3], [1, 2, 3], [1, 4, 2]]
    >>> all([lll(e0, e1) for e0 in test for e1 in test])
    True
    """
    @classmethod
    def add(self, a, b): return (a[0]+b[0], a[1]+b[1])
    @classmethod
    def sub(self, a, b): return (a[0]-b[0], a[1]-b[1])
    @classmethod
    def neg(self, a): return (-a[0], -a[1])
    @classmethod
    def mul(self, a, b): return (a[0]*b, a[1]*b)
    @classmethod
    def dot(self, a, b): return (a[0]*b[0], a[1]*b[1])
    @classmethod
    def calRR(self, a, b): return sqr(a[0]-b[0])+sqr(a[1]-b[1])
    @classmethod
    def norm2(self, a): return sqr(a[0])+sqr(a[1])

    @classmethod
    def getZero(self): return (0.0, 0.0)
    @classmethod
    def getNum(self, num): return (num, num)
    @classmethod
    def getRandom(self): return (rand.random(), rand.random())
    @classmethod
    def genList(self, n): return [(0.0, 0.0)] * n

import random
rand = random.Random(181)

class Vector3(Vector):
    """
    >>> v = Vector3
    >>> a, b = (3, 4, 0), (0, 0, 12);
    >>> v.norm2(v.sub(a, b))
    169
    >>> v.add(a, b) == (3, 4, 12)
    True
    >>> v.neg(a) == (-3, -4, 0)
    True
    >>> v.mul(a, 2) == (6, 8, 0)
    True
    >>> v.dot(a, b) == (0, 0, 0)
    True
    >>> lll = lambda e0, e1: v.norm2(v.sub(e0, e1)) == v.calRR(e0, e1)
    >>> test = [[3, 1, 3], [1, 2, 3], [1, 4, 2]]
    >>> all([lll(e0, e1) for e0 in test for e1 in test])
    True
    """    
    @classmethod
    def add(self, a, b): return (a[0]+b[0], a[1]+b[1], a[2]+b[2])
    @classmethod
    def sub(self, a, b): return (a[0]-b[0], a[1]-b[1], a[2]-b[2])
    @classmethod
    def neg(self, a): return (-a[0], -a[1], -a[2])    
    @classmethod
    def mul(self, a, b): return (a[0]*b, a[1]*b, a[2]*b)
    @classmethod
    def dot(self, a, b): return (a[0]*b[0], a[1]*b[1], a[2]*b[2])
    @classmethod
    def calRR(self, a, b): return sqr(a[0]-b[0])+sqr(a[1]-b[1])+sqr(a[2]-b[2])
    @classmethod
    def norm2(self, a): return sqr(a[0])+sqr(a[1])+sqr(a[2])

    @classmethod
    def getZero(self): return (0.0, 0.0, 0.0)
    @classmethod
    def getNum(self, num): return (num, num, num)
    @classmethod
    def getRandom(self): return (rand.random(), rand.random(), rand.random())
    @classmethod
    def genList(self, n): return [(0.0, 0.0, 0.0)] * n

class Vector4(Vector):
    @classmethod
    def add(self, a, b): return (a[0]+b[0], a[1]+b[1], a[2]+b[2], a[3]+b[3])
    @classmethod
    def sub(self, a, b): return (a[0]-b[0], a[1]-b[1], a[2]-b[2], a[3]-b[3])
    @classmethod
    def neg(self, a): return (-a[0], -a[1], -a[2], -a[3])
    @classmethod
    def mul(self, a, b): return (a[0]*b, a[1]*b, a[2]*b, a[3]*b)
    @classmethod
    def dot(self, a, b): return (a[0]*b[0], a[1]*b[1], a[2]*b[2], a[3]*b[3])
    @classmethod
    def calRR(self, a, b): return sqr(a[0]-b[0])+sqr(a[1]-b[1])+sqr(a[2]-b[2])+sqr(a[3]-b[3])
    @classmethod
    def norm2(self, a): return sqr(a[0])+sqr(a[1])+sqr(a[2])+sqr(a[3])



class Sigma(object):
    """
    >>> s=Sigma()
    >>> print s.add(1,2,3)(4)(5).add(6)
    Sigma: sum is 21, num is 6
    """
    def __init__(self, sig=lambda a,b:a+b, sum=0, n=0):
        self.sig = sig
        self.sum = sum
        self.n = n
        
    def __call__(self, value):
        self.sum = self.sig(self.sum, value)
        self.n += 1
        return self

    def add(self, *value):
        for v in value:
            self.__call__(v)
        return self

    def reset(self):
        self.sum = 0
        self.n = 0
    
    def result(self): return self.sum
    def num(self): return self.n

    def __str__(self): return "Sigma: sum is {}, num is {}".format(self.sum, self.n)
    #def __repr__(self): return "Sigma({},{},{})".format(repr(self.sig), self.sum, self.n)

if __name__ == "__main__":
    import doctest
    doctest.testmod()
